Using propagation for solving complex arithmetic 

constraints 



M.H. van Emden and B. Moa 

Department of Computer Science, 
University of Victoria, Victoria, Canada 
{vanemden, bmoa}@cs . uvic . ca, 
WWW home page: http : / /www . cs . uvic . ca / "vanemden/ 
April 4, 2003; as submitted to CP'03 



Abstract. Solving a system of nonlinear inequalities is an important problem for 
which conventional numerical analysis has no satisfactory method. With a box- 
consistency algorithm one can compute a cover for the solution set to arbitrarily 
close approximation. Because of difficulties in the use of propagation for complex 
arithmetic expressions, box consistency is computed with interval arithmetic. In 
this paper we present theorems that support a simple modification of propagation 
that allows complex arithmetic expressions to be handled efficiently. The version 
of box consistency that is obtained in this way is stronger than when interval 
arithmetic is used. 



1 Introduction 

One of the most important applications of constraint programming is to solve a system 
of numeric inequalities: 

gi(xi , x 2 , ■ ■■ , x m ) < 
g 2 {x x , x 2 , ... , x m ) < 

9k{x\ ,x 2 ,..., x m ) < 

Such systems often appear as conditions in optimization problems. Inequalities are 
regarded as intractable in conventional numerical analysis. The Kuhn-Tucker conditions 
allow these to be converted to equalities. The continuation method |9| is a fairly, but 
not totally, dependable method for solving such equalities. Moreover, it is restricted to 
polynomials. 

An important contribution of constraint programming is the box-consistency method 
1 3 7 1, which improves on the continuation method in several ways. It applies not only 
to polynomials gi,...,gk but to any functions that can be defined by an expression that 
can be evaluated in interval arithmetic. It computes a cover for the set of solutions and 
can make it fit arbitrarily closely. In this way, all solutions are found and are approxi- 
mated as closely as required. The performance of the box-consistency method compares 
favorably with that of the continuation method |7). 
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2 Preliminaries 

In this section we provide background by reviewing some basic concepts. 

2.1 Constraint satisfaction problems 

In a constraint satisfaction problem ( CSP) one has a set of constraints, each of which 
is an instance of a formula. Each of the variables in the formula is associated with 
a domain, which is the set of values that are possible for the variable concerned. A 
solution is a choice of a domain element for each variable that makes all constraints 
true. 

With each type of constraint, there is an associated domain reduction operator, DRO 
for short. This operator may remove from the domains of each of the variables in the 
constraint certain values that do not satisfy the constraint, given that the other variables 
of the constraint are restricted to their associated domains. Any DRO is contracting, 
monotonic, and idempotent. 

When the DROs of the constraints are applied in a "fair" order, the domains con- 
verge to a limit or one of the domains becomes empty. A sequence of DROs activations 
is fair if every one of them occurs an infinite number of times 1 1 01 1 1 . The resulting 
cartesian product of the domains becomes the greatest common fixpoint of the DROs 
1 1 01 1 1 . If one of the domains becomes empty, it follows that no solutions exist within the 
initial domains. This, or any variant that leads to the same result, is called a constraint 
propagation algorithm. 

In practice, we are restricted to the domains that are representable in a computer. As 
there are only a finite number of these, constraint propagation terminates. 

2.2 Constraint propagation 

A Generic Propagation Algorithm (GPA in the sequel) is a fair sequence of DROs. A 
GPA maintains a pool of DROs, called active set, that still need to be applied. No order 
is specified for applying these operators. Even though many variants of GPA exist (see 
Apt 1 1 1 and Bartak 1 2 1), they are all similar to the pseudo-code given in Figure^ Notice 
that the active set A is initialized to contain all constraints. 

put all constraints into the active set A 
while ( A ^ 0) { 

choose a constraint C from the active set A 

apply the DRO associated with C 

if one of the domains has become empty then stop 

add to A all constraints involving variables whose domains have changed, if any 
remove C from A 

} 



Fig. 1. A pseudo-code for GPA. 
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2.3 Interval constraint satisfaction problems 

The constraint programming paradigm is very general. It applies to domains as different 
as booleans, integers, finite symbolic domains, and reals. In this paper we consider 
interval CSPs, which are CSPs where there is only one type and it is equal to the set 1Z 
of real numbers. In such CSPs domains are restricted to intervals, as reviewed below. 

2.4 Intervals 

A floating-point number is any element of F U {— oo, +oo}, where F is a finite set of 
reals. A floating-point interval is a closed connected set of reals, where the bounds, in 
so far as they exist, are floating-point numbers. When we write "interval" without quali- 
fication, we mean a floating-point interval. A canonical interval is a non-empty interval 
that does not properly contain an interval. For every finite non-empty interval X, lb(X) 
and rb(X) denote the left and right bound of X respectively. For an unbounded X, the 
left or right bound is defined as — oo or +oo, or both. Thus, X = [lb(X), rb(X)] is a 
convenient notation for all non-empty intervals, bounded or not. 

If a real x is not a floating-point number, then there is a unique canonical interval 
containing it. Otherwise, there are three. Either way, there is a unique least canonical 
interval containing x, and it is denoted x. 

A box is a cartesian product of intervals. 

2.5 Box consistency 

In 1 7 1, box consistency is computed by a relaxation algorithm implemented in interval 
arithmetic. The algorithm takes as input certain intervals X\ , . . . , X m for the variables 
xi, . . . , x m . It uses each of the functions g\, . . . , in the way that is explained by a 
generic instance that we temporarily call g. We assume that the function g is defined by 
an expression E containing no variables other than x\, . . . , x n . We call the algorithm 
in 1 7 1 a relaxation algorithm because it improves the intervals for the variables one at 
a time while keeping the intervals for all the other variables fixed. This is similar to the 
relaxation algorithms of numerical analysis. 

For simplicity of notation, let us assume that the interval for x\ is to be improved on 
the basis of the fixed intervals X2, ■ ■ ■ , X m for the variables X2, ■ ■ ■ , x m . This is done 
by means of a function gx 2 ,...,X m ( x ) tnat i s defined by evaluating in interval arithmetic 
the expression E with x substituted for x\ and X 2 , ■ ■ . , X m substituted for x%, ... , x m , 
respectively. Thus, gx 2 x m maps a real to an interval. 

To improve the interval X\ = [lb(Xi), rb(Xi)] for x\, suppose that for some a < 
rb(Xi) we have that 

lb(g X2 ,...,x m ([a,rb(X 1 )]))>0. (2) 

In that case the interval for x\ can be improved from X\ to [lb(Xi), a}. 

A bisection algorithm is used to find the least floating-point a for which (0 holds, 
for fixed intervals X2, ■ ■ ■ , X m . A similar bisection is used to improve the lower bound 
of X\ using g and the fixed intervals X2, ■ ■ ■ , X rn . This exhausts what can be done with 
g and X2, ■ ■ ■ , X m . In general, repeating this process with the other arguments and with 
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the other functions among {gi, . . . , gk} causes reductions of X2, ■ ■ ■ , X m and further 
reductions of Xi . 

If one of the intervals becomes empty, it is shown that no solution exists within the 
original intervals X\,... ,X m . Otherwise, the box-consistency algorithm terminates 
with no interval reduction possible. Let us call the resulting state of the domains func- 
tional box consistency, to distinguish it from the type of box consistency described 
below. 

As the criterion for a being an improved upper bound for X\ is (0 with the left-hand 
side evaluated in interval arithmetic, this box consistency algorithm can be improved 
by means of interval constraints, as was pointed out in 111 II . Here it was proposed that 
instead of such interval arithmetic evaluation, one applies propagation to the interval 
CSP containing as constraints 

g{x\, . . .,x m ) < 

x\ > a (3) 
x 2 e x 2 , . . . , X rn e X rn 

The result is called relational box consistency. In II II it is shown that the resulting 
intervals are contained in those obtained by functional box consistency. 

To apply propagation, one needs to decompose the arbitrarily complex expression 
for g into multiple primitive arithmetic constraints (x + y — z, x x y = z, as well as 
those involving trigonometric or logarithmic constraints), as explained in section [5] so 
that the corresponding efficient DROs can be applied. In this way the structure of E is 
lost. As a result, GPA will activate DROs that have no effect, even though, on eventual 
termination of GPA, the result is the correct one: the unique and consistent state, or 
failure. 

This problem was addressed in 1411 II . The remedy described there is to create a 
tree data structure for E and perform propagation based on the tree structure. Such a 
structured propagation algorithm avoids superfluous activations of DROs by following 
the tree from the bottom to the top and then from the top to the bottom and to repeat 
these two traversals as a cycle. 

In this paper we show that by a simple modification of the GPA one can get the effect 
of an optimized version of a structured propagation algorithm, yet without maintaining 
a tree data structure. The downward part of this algorithm is enhanced by the absence of 
multiple occurrences of variables. In interval arithmetic this would be an unacceptable 
limitation. However constraint programming allows us to translate multiple occurrences 
to equality constraints. This may be a welcome improvement compared to the usual way 
of representing a system such as in Equation^ 

In section|3] we describe our translation of a system such as in Equation[2to a CSP 
consisting of primitive constraints. 

3 Generating a CSP from a system of nonlinear numerical 
inequalities 

The functions gi, . . . , gk in Equation^are defined by expressions that can be evaluated 
in interval arithmetic and hence can be translated to CSPs containing only primitive 
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constraints. Any of these expressions may have multiple occurrences of some of the 
variables. As there are certain advantages in avoiding multiple occurrences of variables 
in the same expression, we rewrite without loss of generality the system in Equation^ 
to the system shown in Figure|2] 



gi(xi , X2 , ■ ■ ■ , x m ) < 

02 (ail , X2 , ... , %m) < 



g k (xi , x 2 , • ■ • , x m ) < 
al/Eq (vi s i , Ui,2 , • ■ • , Ui,m) 

allEq (v Pt i , v Pt 2 , • ■ • , v P ,n p ) 

Fig. 2. A system of non-linear inequalities without multiple occurrences. The variables 
{#i, . . . , x m } are partitioned into equivalence classes V\, . . . , V p where Vj is a subset 
{v jA , . . . ,v jt7lj } of {xx, ■ ■ -,x m }, for j 6 {1, . . . ,p}. 



In Figure |2] the expressions for the functions <?i , • • • , <?fc have no multiple occur- 
rences of variables. As a result, they have up to m rather than n variables, where m > n. 
This special form is obtained by associating with each of the variables Xj in Equation^ 
an equivalence class of the variables in Figure |2] In each expression each occurrence 
of a variable is replaced by a different element of the corresponding equivalence class. 
This can be done by making each equivalence class as large as the largest number of 
multiple occurrences. The predicate allEq is true if and only if all its real-valued argu- 
ments are equal. 

An advantage of this translation is that evaluation in interval arithmetic of each 
expression gives the best possible result, namely the range of the function values. At 
the same time, the allEq constraint is easy to enforce by making all intervals of the 
variables in the constraint equal to their common intersection. This takes information 
into account from all k expressions. If the system in its original form[2 with multiple 
occurrences, would be translated to a CSP, then only multiple occurrences in a single 
expression would be exploited at one time. 

The translation of a complex arithmetic expression to a CSP containing primitive 
arithmetic constraints is an obvious variant of the procedure that has been familiar since 
FORTRAN compilers parsed a complex arithmetic expression and generated code from 
the parse. The CSP variant of this procedure is at least as old as BNR Prolog 1 5 1, which 
was first implemented in the late 1980s. For a formal description of the translation we 
refer to [ 1 1 1. We give a brief informal description here. 

For the purpose of the translation, one regards an expression as a tree with operators 
as internal nodes and constants or variables as external nodes. We associate with each 
internal node a unique variable. Each internal node now generates a primitive constraint. 
For example, if an internal node has "/" as an operator, a; as a variable, and y and z 
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as variables associated with left and right child nodes, respectively, then the ternary 
constraint x x z = y is generated. 

In this way, each expression generates a primitive constraint for each internal node. 
Let v be the variable associated with the root. As this represents the value of the entire 
expression, the primitive constraint v < is generated as well. Because of the absence 
of multiple occurrences, the partial CSP generated is an acyclic CSP for which local 
consistency implies global consistency |8|. 

Finally, as the allEq constraints in Figure|2]do not contain expressions, they need no 
translation to more primitive constraints: they have the obvious, optimal, and efficient 
domain reduction operator described earlier. 

This completes our description of how to translate a system as in Equation ^ to a 
CSP. To solve it by propagation, we must first consider propagation in a CSP consisting 
only of primitive constraints generated from the tree of a single expression together with 
the inequality constraint involving the root variable. This we do in the next section. 

Before proceeding thus, we point out that translating the system <|3J to a CSP in the 
way just described enhances the opportunities for parallelism in propagation beyond 
those already present in the system Q. To investigate these would take us beyond the 
limits of this paper. However, it will be useful here to highlight the structure that gives 
rise to these opportunities by means of a hardware metaphor. 

In the first place it is important to note that the sets of constraints arising from the 
same expression form a cluster for the purposes of propagation. For example, with the 
exception of the root, the internal variables only have unique occurrences. As a result, 
when the DRO of a constraint is activated, it usually causes constraints generated by 
the same expression to be added to the active set. However, the external variables may 
occur in all of these clusters. 

For the purposes of a parallel algorithm it is useful to imagine the clusters arising 
from each of the expression as hardware "cards", each connected to a "bus", where 
the lines of the bus represent the external variables in common to several expressions. 
Whenever two external variables belong to the same equivalence class, they are con- 
nected by a "jumper" in the hardware model. 

The hardware architecture suggests a parallel process for each card that asynchronously 
executes DROs of constraints only involving internal variables. The processes synchro- 
nize when they access one of the bus variables. The DROs of allEq constraints can also 
be executed by a parallel process dedicated to each. 

4 Modifying propagation for evaluating an expression 

We first show that, regardless of efficiency, GPA can be used to evaluate a single ex- 
pression. Suppose that we have an expression E in variables x\, . . . , x n . Let y be the 
variable at the root of the tree representing the value of E. Let C be the CSP generated 
by E as described before. 

Theorem 1. Suppose the domains ofxt, . . . the intervals X\ , . . . . X n . Suppose 

the domains of y and the other internal variables are [— oo, +oo]. Applying the GPA to 
C results in the domain of y being the same interval as the one obtained by evaluating 
E in interval arithmetic with Xi,..., X n substituted for Xx, . . . , x n , respectively. 
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Proof. According to I1QI1I . every fair sequence of DROs in GPA converges to the 
same limit for the domains of the variables. There is a finite sequence s of DROs that 
mimics the evaluation of E in interval arithmetic. At the end of this, y has the value 
computed by interval arithmetic and GPA terminates. 

Theoremnjis useful in showing that, in the absence of information about the value of 
the expression, propagation does the equivalent of interval arithmetic. But the GPA does 
it in a wasteful way. GPA does not specify the order of applying the DROs other than 
that their sequence should be a fair one. In a typical random fair sequence, many DROs 
will not have any effect. This inefficient behavior is the motivation for our modifications 
to GPA presented here. 

Theorem 2. Suppose we modify GPA so that the active set is initialized to contain 
instead of all constraints only those containing at most one internal variable. Suppose 
also that the active set is a queue in which the constraints are initially ordered according 
to the level they occupy in the expression tree, with those that are further away from the 
root placed nearer to the front of the queue. Then GPA terminates after activating the 
DRO of every constraint at most once. On termination, y has as domain the value that 
the expression has in interval arithmetic. 

Thus we see that GPA has exactly the right behavior for the evaluation of an expres- 
sion if only we initialize the active set with the right selection of constraints. We call 
this propagation with selective initialization (PSI). In the sequel, the constraints that 
have at most one internal variable are referred to as peripheral constraints. 

5 Selective initialization for obtaining box consistency 

Propagation, whether modified or not, obtains results that are at least as strong, and 
typically stronger than, box consistency as described in 1 7 1 . This can already be demon- 
strated when considering a single expression E in variables x±, . . . , x m with intervals 
X\, . . . , X m as domains. A step towards box consistency is to reduce the interval for 
each of the variables separately. To simplify notation we do this for x\, keeping the 
interval domains X2, ■ ■ ■ , X m for x-x, ... , x m fixed. Suppose that for some a < lb(Xi) 
we have that lb(gi([a, rb(Xi)], X2, ■ ■ ■ , X m )) > 0. In that case the interval for X\ can 
be improved from X\ to \lb(X\), a]. 

Suppose that instead of such an interval arithmetic evaluation, one applies propaga- 
tion to the interval constraint system containing as constraints x\ > a, X2 £ X2, ■ ■ ■ , x m 6 
X m , as well as all the ones obtained by translating the expression tree of E to primitive 
constraints. 

Theorem 3. Suppose that lb(E([a, rb(X{j\, X 2 , . . . ,x m £ X m ) > 0. When GPA is 
applied to this CSP, failure results. 

Proof. Consider any fair sequence that starts with a segment s mimicking the inter- 
val arithmetic evaluation of E. At the end of s, the interval for y has a positive lower 
bound, by the assumption. The fair sequence can be continued by applying the DRO for 
y < 0. This yields failure. 
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This proves the theorem. If we use GPA in this way for box consistency instead 
of interval arithmetic, we never obtain a worse result and we typically obtain a better 
result. 

However, unmodified GPA will obtain the better result in an inefficient way. 
Our next result is a modification of GPA that obtains the better result in a more 
efficient way. 

Suppose we have a CSP S generated by an expression E. Let y be the variable at the 
root of the tree representing E. Suppose we apply GPA to S. After the termination of 
GPA, we have the following theorem. 

Theorem 4. Suppose the domain for y is changed to a proper subset of it. Suppose the 
unique constraint containing y is placed in the active set as only element. Then GPA 
terminates with the same result as when the active set would have been initialized to 
contain all constraints. 

Proof. Let s be any fair sequence of DROs. If s starts with the unique constraint con- 
taining y, then the theorem follows. If s starts with a different constraint c, then applying 
the DRO associated with c does not affect any domain (DROs are idempotent and the 
DRO of c is already in its fixpoint). Thus, removing c from s does not affect the fix- 
point of s. We keep removing the first element of s until we reach the unique constraint 
containing y. The new s' formed starts with the unique constraint containing y and has 
the same fixpoint as s. Thus, Theorem|4]is proved. 

This theorem shows that a complex constraint of the form E(x%, . . . , x n ) < can 
be made relationally box consistent by using a modification of GPA that is more efficient 
without affecting the quality of the result. 

The result obtained from Theorem|4]is usually better than backward evaluation |4 1. 
The following example illustrates that. 

Let us consider the following system. 

xi/x 2 < 
xx e [-1,1] x 2 e[-l,l] 

Using the decomposition described in section|3] we generate the following CSP. 

y < o 
x\/x 2 = y 
x\ e [— l, 1], X2 e [— l, 1], y e [-co, +oo] 

Applying the GPA described in Theorem |4] gives the fixpoint [—1, 1] x [—1, 1] x 
[—oo, — 1]. The result obtained from interval arithmetic backward evaluation is [—1, 1] x 
[-1,1] x [-oo,0]. 

Pseudo-code for the PSI algorithm is given in Figure|5] 

In certain situations (fully described in [ 12 1), Theorem|5] stated below, can be used 
instead of Theorem|4] 

In addition to the assumptions of Theorem^] we suppose that the expression E has 
no multiple occurrences of any variable. 



9 



put only peripheral constraints into the active set A 
while ( A ^ 0) { 

choose a constraint C from the active set A 

apply the DRO associated with C 

if one of the domains is empty then stop 

add to A all constraints involving variables whose domains have changed, if any 
remove C from A 

} 

Fig. 3. Pseudo-code for Propagation by Selective Initialization. 

Theorem 5. Suppose the domain for y is changed to a proper subset of it. Suppose the 
unique constraint containing y is placed in the active set as only element. Then GPA 
terminates after having activated the DRO of each constraint at most once. 

Moreover, the result is the same as when the active set would have been initialized 
to contain all constraints. 

This theorem is based on the fact that once the fixpoint is obtained, reducing the 
domain of one variable may cause the other domains to be reduced but not itself. As 
shown in 1 12], this is true when DROs have certain properties. For example, one could 
allow domains to be the union of disjoint intervals, as in the systems of Hyvonen or 
Havens 18 161 . But when DROs are those described in this paper, reducing the domain 
of a variable can affect its own domain as shown in the example above. Even though 
[—1, 1] x [—1, 1] x [—oo, +oo] is a fixpoint of the CSP 

xi/%2 = y 
xi e [— 1, 1], x 2 e [-1, l], y e [— oo, +oo] 

reducing the domain of y to [— oo, 0] leads to a further reduced domain of y that is equal 
to [— oo, —1]. 

6 Conclusion 

In this paper, we presented a slight modification to propagation to get the same or better 
results than structured propagation. Such a modification is used in a box-consistency al- 
gorithm to solve systems of non-linear inequalities but it can also be used when solving 
other CSPs obtained by translating complex expressions. In fact, most CSPs of practical 
importance seem to be derived from complex expressions. 
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